Backflow Correlations for the Electron Gas and Metallic Hydrogen 
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We justify and evaluate backflow-threebody wavefunctions for a two component system of elec- 
trons and protons. Based on the generalized Feynman-Kacs formula, many-body perturbation 
theory, and band structure calculations, we analyze the use and the analytical form of the backflow 
function from different points of view. The resulting wavefunctions are used in Variational and 
Diffusion Monte Carlo calculations of the electron gas and of solid and liquid metallic hydrogen. For 
the electron gas, the purely analytic backffow and three-body form gives lower energies than those 
of previous calculations. For bcc hydrogen, analytical and optimized backflow-threebody wavefunc- 
tions lead to energies nearly as low as those from using LDA orbitals in the trial wavefunction. 
However, compared to wavefunctions constructed from density functional solutions, backflow wave- 
functions have the advantage of only few parameters to estimate, the ability to include easily and 
accurately electron-electron correlations, and that they can be directly generalized from the crystal 
to a disordered liquid of protons. 
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I. INTRODUCTION 

This paper concerns the form of the ground state wave- 
function of metallic hydrogen at high enough density so 
that all the hydrogen molecules are dissociated and the 
electrons are delocalized. Neglecting possible quantum 
effects on the protonic motion, the many-body wavefunc- 
tion can be regarded as the ground state of an electron 
gas under the influence of an external potential due to the 
actual positions of the protons. Quantum Monte Carlo 
(QMC) techniques are currently one of the most power- 
ful methods to calculate accurately the properties of such 
a many-body quantum systemi. However, since ground 
state QMC is based on trial wavefunctions, QMC typi- 
cally demands compact and accurate descriptions of the 
ground state wavefunction. In this paper we review dif- 
ferent approaches to obtain and improve trial wavefunc- 
tions, compare the qualities of the resulting many-body 
wavefunctions with previous QMC calculations for the 
electron gas and metallic crystal hydrogen, and present 
first results using these wavefunctions for liquid metallic 
hydrogen. 

Most of the work within QMC has been done using a 
pair product (PP) (or Slater-Jastrow) wavefunction: a 
Slater determinant of single electron spin orbits times a 
product of pair electron (Jastrow) factors. Notwithstand- 
ing certain deficiencies such as a lack of direct spin cou- 
pling, this wavefunction has proven to be quite accurate, 
in particular within fixed-node Diffusion Monte CarloA 
(DMC). The first calculation on many-body hydrogen^, 
used an even simpler form of this wavefunction; the sin- 
gle electron orbits were taken to be free electron plane 
waves. We refer to this as the SJ-PW trial function. 
Later, Natol:'^-^ found that determinants using these or- 
bitals are inaccurate by 0.05eV/atom within the fixed- 
node DMC calculations at the density corresponding to 



the transition between molecular and metallic hydrogen 
{vs = 1.31). Hence, more accurate orbitals, computed 
from either density functional (LDA) or Hartree-Fock 
(HF) calculations, are required. Because these orbitals 
are calculated assuming fixed ionic positions, inclusion 
of ionic motions, such as those from the zero point mo- 
tion of the ions in the crystal, is difficult. 

Recently, there have been new attemptE'"^'^ to calculate 
properties of disordered systems such as liquid hydrogen 
within QMC. In the Coupled Electron Ion Monte Carlo 
(CEIMC) methoc^ the protons are moved based on the 
results of a QMC calculation of the electronic energy. 
This approach requires accurate trial functions that can 
be obtained quickly as the ionic positions are changed; 
methods involving the solution of mean field equations 
such as LDA and HF, or even optimizing a parameterized 
trial function, can greatly slow down the overall perfor- 
mance of the CEIMC simulation'^. Further, combining 
the orbitals obtained from LDA or HF with a pair cor- 
relation (Jastrow) factor to improve the accuracy is not 
straightforward; substantial modification of the orbitals 
might be necessary requiring a reoptimization of the or- 
bitals and correlation factor, in principle, at each ionic 
position^. This optimization step creates a bottleneck 
to coupling the QMC calculations with the ionic Monte 
Carlo. 

One could consider obtaining the trial wavefunction 
from other variational approaches like Fermi-hypernetted 
(FHNC) chain or correlated basis functions (CBF) 
methods^ which would not have the problems of opti- 
mization. However, in these approaches based on explicit 
integration, one is in general limited in the form of the 
trial function by the ease performing the integration, and 
these are typically much more time-consuming than LDA 
calculations. 

One of the biggest advantages of the QMC approach 
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is that one can use an arbitrary wavefunction without 
changing the algorithm in an essential way. Fast algo- 
rithms will result if one can find concise and accurate 
forms. In this paper, instead of using one-body orbitals 
from mean field theory or integral equations, we propose 
to use trial functions which depend explicitly and contin- 
uously on the ionic variables. Such wavefunctions do not 
have to be reoptimized for movements of the ions, are 
easy to implement, and accurate for disordered systems. 
Calculation of ionic forces is simplified since the deriva- 
tive of the trial function with respect to ionic configura- 
tions is a straightforward application of the chain rule. 
These trial functions are a generalization of the back- 
flow three-body wavefunctions used very successfully in 
highly correlated homogeneous quantum liquids: liquid 
■^He and the electron gas. There, backflow trial functions 
show much improvement over the pair product getting 
approximately 75% of the energy missing at the PP level 
and even more when done with the fixed-node method. 

Backflow wavefunctions were developed by Feynman 
and CohenS for a single ''He impurity in liquid ''He when 
it was recognized that without backflow, the mass of the 
impurity was equal to the bare mass. Pandharipande and 
Itoh'^" showed that the backflow arises from the momen- 
tum dependence of the correlation between the impurity 
and the liquid. The backflow wavefunction was then ex- 
tended to bulk liquid ^Heiiti^ using an integral equation 
method to evaluate expectation values. The first use of 
backflow in QMC was by Lee et alm^ and othersi^^ with 
calculations on liquid ^He. Moroni et ali^ further op- 
timized the trial function within liquid ^He. Kwon et 
ali^L^ used backflow functions for the electron gas in 
both 2D and 3D, obtaining significantly lower energies 
and improved excitation energies. Vitiello et al»^ dis- 
cuss an equivalence of backflow and spin-dependent cor- 
relations, an aspect we will not further consider in this 
paper. 

Using different approaches, we generalize the back- 
flow three-body wavefunction to a two component sys- 
tem of electrons and protons and derive approximate 
expressions for the correlated trial function. We first 
present an argument based on the generalized Feynman- 
Kacs formula which shows that backflow is the next order 
improvement beyond the pair product (PP) wave func- 
tion. Using perturbation theory, we then discuss gen- 
eral features of the backflow functions and obtain ex- 
plicit expressions for the homogeneous electron gas and 
for the electron-proton plasma. A similar analysis us- 
ing the Bohm-Pines method has been recently performed 
by Gaudoin et al.^", however, without going beyond the 
Slater- Jastrow wavefunction. Studying the problem of 
a single electron in the potential generated by a simple 
cubic lattice of protons, we show that the exact one elec- 
tron wavefunction can be approximately rewritten by a 
backflow function. Finally, we optimize numerically sim- 
ple functional forms for the backflow functions in the 
full many-body problem by variational Monte Carlo. We 
compare the quality of the wavefunctions stemming from 



these different approaches for the electron gas and for liq- 
uid and crystal hydrogen at the level of variational and 
diffusion Monte Carlo. 

In the following we consider the non-relativistic Hamil- 
tonian of TV protons and N electrons: 

where Ai = f? ji^nii), i — 1,...,2N and and e, 
are the electron or proton mass and charge. The Fermi 
wavevector is kp- Numerical results are given in atomic 
units where Ag = 1/2 and Ap = for classical pro- 
tons, 1 6,1 = rrii = 1. The electron density n = N/V 
is quoted in terms = a/og, where a = {4Tm/3)~^^^ 
and Go — h^/mf.e^ is the Bohr radius. Energies of the 
QMC calculations are given in Rydbergs per electron. 

II. THE FEYNMAN-KACS APPROACH TO 
IMPROVING THE WAVEFUNCTION 

The Feynman-Kacs formula expresses the exact wave- 
function in terms of average over Brownian paths. We 
now review how it can be generalized to random walks 
with "drift". 

We define the "importance-sampled" Green's function 
as Gt — ip exp{—tH)ip^^ in operator notation where -0 is 
an unsymmetrical trial function. Gt acting on a function 
has the effect of enhancing the component of lower energy 
states. Then the lowest energy (exact) fermi wavefunc- 
tion (j>F{R) is given by: 

(j)F{R) oc A^p{R) \im I dR'{R'\Gt\R) (2) 

t^oo J 

assuming only that the trial function has a non-zero over- 
lap with (f)F and that (ftp non-degenerate. Here A is a 
projection operator for fermion symmetry defined as 

M{R)^^X.(-^^''f(P^^ (3) 
p 

and R = {ri,r2,...} is a point in configuration space. 
The electron spin is treated by restricting the permuta- 
tion in Eq. (|3J) to be exclusively within spin up or spin 
down electrons. 

Following the derivation in Diffusion Monte CarloSi, 
the Green's function can be split into diffusion, drift and 
branching processes. To show this, the master equation 
for the Green's function is written: 

dG 

= ^iJV'-^G = A,V,(V,-H2VanV')+£;(i?)]G. 

i 

(4) 

The local energy E, defined as E{R) = ip'^Hip, is the 
residual error of the trial function, and becomes a con- 
stant function as ijj approaches an exact eigenvalue. Trot- 
ter's formula applies to the above master equation, allow- 
ing us to split up the evolution into the first two terms 
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describing a stochastic process, and the final term which 
is a branching or "weighting" process. Thus we have the 
generahzed Feynman-Kacs formula : 



bpiR) oc A4iR) ^exp 



dtE{R{t)) 



(5) 



where the brackets imply averaging over all drifting ran- 
dom walks R{t) beginning at a point R. The above rela- 
tion is exact for any real trial function. For trial functions 
having an imaginary component of V Ini/;, the formalism 
goes through, however, the Green's function is no longer 
real and positive and therefore cannot be treated as a 
probability. Other methods are more appropriate. For 
the moment we will ignore this case. 

To make further analytical progress, we take the aver- 
age into the exponent. For any stochastic process, one 
can write the average of the exponent as the exponential 
of the cumulant expansion, the first two terms of which 



(j>F{R) ^ A'rP{R)eM-{{E)) + [l/mSE'')) 



(6) 



The double brackets are defined as {{E)) = 
(/p dtE{R{t))) with walks R{t) generated from the drift 
and diffusion starting at a point R. We truncate the 
cumulant expansion after the first term. We then have 
an approximate method of improving the trial function. 



(7) 



with the subscript indicating that the drift is given by 
V Ini/i^"). If we split the log of the trial function into its 
real and imaginary parts V'*'"' = exp(— C/*^") -l-iS'^"') with 
U and S real, we are led to the following equations for a 
single iteration: 

jj(n+l) ^ jj(n) 

+ (( V +^A,[V2c/(")-(V,C/("))2 + (V,5("))2))„ 



(8) 



+ 1) 



(9) 

Here V{R) is the total potential energy. 

Specializing to the case of a fermi liquid, we take as an 
initial wavefunction [7*^°^ = and S*'"' = X^jl^i ' 
singly occupied free particle states. (The usual spin func- 
tions are assumed but not explicitly written.) Note that 
this function is an unsymmetrical trial function, with 
a non-zero overlap with a fermion state as long as all 
the ki's are distinct. When the wavefunction is antisym- 
metrized, one gets a determinant of plane waves. How- 
ever, the antisymmetrization will be done only once, after 
the trial function has gone through several iterations of 
Eq. 0. This way simplify the procedure, since the local 
energy of the unsymmctric trial function is much simpler 



than that of an antisymmetric trial function. Note that 
in Eq. (5) both the antisymmetrization and the averag- 
ing are linear operators and so can be interchanged. 

After the first iteration, the wavefunction will have the 
form: 



[/(i^ = ((V(i?)))o ^ U{R) 



(10) 



(11) 



In the above equation and the following discussion we 
drop, without mention, constant normalization terms. If 
V{R) — '^i^jv{rij) is a pair potential, with a Fourier 
transform Vk, the averaging can be carried out analyti- 
cally with a result that C/'^-' will also be a pair poten- 
tial and will have a Fourier transform given by Vk/ (Afc^) 
where A — Xi + Xj. For a Coulomb potential the real- 
space correlation (Jastrow) function will then have the 
form: Mc(r) — —e^r/2X. 

Hence, the form of the first order wavefunction is 
of the SJ-PW or Slater- Jastrow form, with free parti- 
cle orbitals. The pair term will be denoted by U with 
U — ^i^jU(rij). Typically the form of u is derived 
from a variational principle, chosen so such that either 
the total energy or variance is minimized. This will, of 
course, give a lower energy than the cumulant form de- 
rived above. The above derivation does give the correct 
cusp condition (the limit of u at large k or small r) . How- 
ever, it does not give the long wavelength limit correctly 
because of the neglect of the higher cumulants. GaskellS^ 
proposed an analytic form based on the Random Phase 
Approximation (RPA) without any parameters. It was 
found^'^ for the homogeneous electron gas that the RPA 
form does, as well as, or better than simple assumed 
forms with parameters. Figure ^ shows a comparison 
of these correlation functions. 

Note that the cumulant approximation will not exist 
if the Fourier transform of the potential does not ex- 
ist. Two examples of such potentials are the hard sphere 
and Lennard- Jones interactions. However, for the short 
range part of a soft potential which does have a Fourier 
transform such as the Yukawa potential, the cumulant 
approximation works quite well (see remarks concerning 
the situation at finite temperature in Ref.^'*). 

We now perform the next iteration of this procedure. 
To minimize the fluctuations in the local energy so that 
the cumulant approximation will be more accurate, we 
assume that first order wavefunction has been optimized 
but it still has a pair product form. Using Eq. lO, 
neglecting constants and combining pair terms together, 
we get in second order a function of the form: 

f/(2)=^(i)_((^A.(V.C/W)2))i (12) 



and 



5(2) = ((^k,. (r,-2A,V,{/)))i 



(13) 
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tion term W. The backflow potential is 
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FIG. 1: The electron-proton Jastrow 
{aoG)~^ from band calculations of solid cubic hydrogen at 
Vs = 1.31 (squares), Eq. (53). The rightmost square is 
the first reciprocal lattice vector. This is compared with the 
RPA (Gaskell) form (solid line) Eq. (59) and cumulant form 
(dotted line) 1/G* and the improved analytic form (Eq. 71) 
(dashed line). 



Here U^^^ includes additional pair terms. 

At second order, we cannot perform the averaging an- 
alytically, since it involves drift under the influence of the 
first order wavefunction: U'^^^ We make the assumption 
that the averaging will not change the functional form of 
the quantity being averaged but only smoothes out the 
individual functions. That is, our ansatz for the iterated 
wavefunction is: 



and 



[/(2) ^[/(i?)_^(V,Ty)' 



5(2) ^^k,.(r,-V,F) 



(14) 



(15) 



where U, W and Y are three different pair "poten- 
tials" to be optimized. In the following, we have adopted 
the convention that pair functions have the same sign 
as Vij{r), so that, for example, a repulsive v leads to a 
repulsive w and y. 

The two new functions appearing at second order are 
the backflow function Y and the three-body or polariza- 



(16) 



where y{r) is a spherically symmetric function and the 
sum extends over all pairs of particles, including both 
electrons and protons. The backflow displacement is de- 
fined as the gradient of the backflow potential with re- 
spect to a particle coordinate: 



Ar,; 



where 



1 dyjr) 
r dr 



(17) 



(18) 



corresponds to the definition in previous work for homo- 
geneous systemsi2ii2i. 

With this ansatz, the antisymmetized trial function is 
a determinant composed of "quasi-particle" coordinates: 



det [expikj(ri + Ar^)] e 



-u{R)+{vwy 



(19) 



Recall that in the fixed-node or fixed-phase diffusion 
Monte Carlo method, one obtains the exact energy sub- 
ject to the imposed constrain'feiiS^. The assumed node or 
phase limits the ultimate accuracy for fermion systems. 
Since the correction to the real part, the three-body term, 
is already symmetric, it is the backflow which is respon- 
sible for the change of node or phase of the trial function 
and is, in that sense, more important than the Jastrow 
and polarization part. 

In the above derivation we have neglected any effects 
of a complex drift velocity. However, as already shown in 
Ortiz and Ceperley2&, a complex drift velocity does not 
affect the corrections to the wavefunction to the order 
we have considered; the Eqs (PJ are valid to improve the 
wavefunction. 

Now we consider the long range properties of the pair 
functions appearing. In periodic boundaries (or "super- 
cells") we need to perform Ewald summations of the 
functions V, U,W^Y. This is most convenient in Fourier 
space. We define the Fourier transform of a radial func- 
tion as: 



(20) 



Using the Poisson sum formula, the "potential" of the i*'* 
particle in periodic boundary conditions is: 



(21) 



where V is the volume of the supercell. For example, to 
find the backflow displacement, Ea. H17|l . we simply take 
the gradient of the pair function: 



Ari 



V 



(22) 
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where k ranges over the reciprocal lattice vectors of the 
supercell. 

The three-body potential, W is defined analogously in 
terms of a pair polarization w{r). This function is related 
to that used in previous QMC worl4i^iiSiiLi& by: 



1 dw{r) 
r dr 



(23) 



The overall sign of w is not important because only its 
square appears in the trial function, but the relative sign 
of the electron-electron to the electron-proton interaction 
is significant. 

One of the simple ways of deriving conditions on the 
backflow function is to look at the action of the Hamilto- 
nian on the wavefunction, the local energy, and to mini- 
mize the fluctuations of the local energy. Here we focus 
on the imaginary part of the local energy and consider a 
single electron with phase S" = q • (r — VY) . Setting to 
zero the imaginary part of the local energy we obtain: 



VV^yW + 2Vu(r) - 2Vu(r)VV?/(r) = 0. 



(24) 



Neglecting the last term, since it is higher order in the 
interaction, we obtain: V^?/(r) = —2u{r). This has a so- 
lution in Fourier space: Hk — 2uk/k'^ ■ (Because we want a 
solution which is smoother than u(r) at r = 0, we neglect 
a term proportional to .) We get the same smoothing 
{k~^) that we observed at first order for the pair function. 
Shown in figure |21 are the u{r) and r]{r) function coming 
from this approach. Note that this approach is based 
on a single electron description and therefore does not 
correctly describe the long wavelength (large r) behavior 
where the collective motion dominates. 

To obtain a simple form for the three-body potential, 
we note that the averages used in the definition of Y are 
similar to those for W, see Ea. (|12ll5|) . Hence an estimate 
of the polarization potential is 



W 



-VXkpY 



(25) 



where we have approximated (((Vit/)^)) « {{V ,U))'^ / t , 
averaged over a "typical" time t sa (Xkp)^^. This relates 
the three-body contribution to the backflow potential. 

The GFK approach is good for suggesting corrections, 
but there are serious problems in using it to find a good 
backflow function since the averaging is difficult to carry 
out, the linear cumulant approximation may be inade- 
quate, and the long-time effects of the imaginary drift 
are being ignored. If one cannot analytically perform the 
averaging, one does not know what time to multiply the 
local energy by to get a wavefunction, nor the relative 
corrections at large versus small distances. We now dis- 
cuss several other approaches which allow us to directly 
evaluate the Jastrow, threebody and backflow functions 
and give more insight into their form. 
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FIG. 2: The u^P{r) using the RPA (Gaskell) form (dotted 
line) and ri{r) (solid line) from smoothing it with for 
the e-p correlation at rs = 1.31, both computed for an infinite 
system. Note that in this approximation they both tend to 
the same limit at large r. 



III. PERTURBATION THEORY/ ANALYTIC 
METHODS 



In this section we follow another approach to obtain- 
ing improved estimates of the many-body wavefunction. 
Many-body perturbation theory is a well studied ap- 
proach to understanding the effects of weak correlation. 
Encouraged by the use of the RPA^^ which gave an ex- 
cellent analytic two-body correlation function, we will 
extend this wavefunction by perturbative expressions for 
the Jastrow, backflow and threebody potentials for the 
electron gas and for metallic hydrogen. Rather than per- 
forming a systematic low or high density expansion to de- 
rive analytical expressions for the variational wavefunc- 
tion of the electron gas or metallic hydrogen, we concen- 
trate on improving this correlation factor. The collective 
coordinate formulation of Bohm and Pines^^'^^ allows us 
to use Slater- Jastrow wavefunctions as zero*''-order start- 
ing point. We obtain improved potentials for the homoge- 
neous electron gas and metallic hydrogen, which compare 
very well with numerically optimized forms. 

Even if perturbation theory assumes a weak coupling 
(or high density) expansion, we expect the derived prop- 
erties to be qualitatively valid as long as the correspond- 
ing perturbation expansion remains regular, e. g. until 
there is a phase transition to an insulating phase. 
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A. Single particle perturbation theory 

Consider a single electron interacting with an arbitrary 
external potential v(r) with Fourier transform w(k). To 
avoid the problems arising from the long-range behavior 
of the Coulomb interaction, we restrict the analysis to a 
potential with a Fourier transform which is finite at the 
origin, |'L'(0)| < oo, e.g. a screened Coulomb potential. 
We use the continuum notation in this section (-^ <-> 

/ (1^). The solution of the Schrodinger equation (/'fc(r) 
of a particle with wavevector k. 



can be written as 



Cfc(p) = (27r)3^(k-p) 



47r/(k,p) 

fc2 — + iS 



(26) 



(27) 



where the off-shell scattering amplitudes, /(k, p), are 
given by the integral equation 



4^A/(k,p)= J ^«(p_k')cfe(k' 



(28) 



Using the Born approximation we can write down the 
wavefunction to first order in v 



Mr) ^ 0r(r)-H0i^Hr) 



1 



v{p) 



(27r)3 p • (p + 2k) 



(29) 



If we expand the solution around k = and assume the 
change in the wavefunction is small, we can write it in 
the pair-product and backflow form, Ea. H19|) . We obtain 
for the pair potential 



u{y) 



(27r)3 

and for the backflow potential 



y(r) 



(27r)3 



Ap2 



Ap4 ■ 



(30) 



(31) 



Note that the small p part of the integral is usually cut 
off by the finite size of the box. In addition, it is ri{r), 
the derivative of y{r) (see Eq. (18)) which enters in the 
trial function. 

Although the first order approximation is only reliable 
in the case of a weak potential, it becomes correct in the 
high momentum region and hence, gives the correct cusp 
conditions. The derived form is identical to that obtained 
from the Feynman-Kacs formula in the previous section. 
For an arbitrary weak potential, we further get the long 
range behavior, u cx v{0)/r and ry oc l/r3 for r — > oo, 
provided the potential has a finite range {v{q) — v{0) cx 
for q —>■ 0) and there is no other singularity in the 
integrand. 



To find an approximate form for the three-body func- 
tion Vl^(r) we must go to higher order in the interaction, 
but only at k = 0. Using Eqs. H27I28|I . we can write 
down the second order corrections in v to the wavefunc- 



tion, 



at fc = 0: 



dPq (Pp w((j)e*'i '' v{p)t 
(27r)3 (27r)3 (q + p)^ p'' 



(32) 



(Pq (fip w(g)e* 
(2^(2^^^ 



1 - 



2q- p\ w(p)e*P ' 



Ap2 



This is almost in the form of the three-body correlation 
obtained with the FK approach: (Vw)^. Note, however, 
that Eq. (33) is unsymmetrical in q and p so that in r- 
space it will be written as: (Wwu) ■ iSwy) with ^^(r) ~ 
u{r) and Wy{r) ~ y(r). Therefore the polarization term 
is not a square but a product of the gradients of two 
different functions. (In second order one will also find a 
contribution oc [u(r)]2 to the pair term.) 

The perturbative expressions (|30|) and (|31|l are based 
on the Born approximation for scattering between free 
states. However, an attractive potential as the electron- 
proton (effective) interaction might also lead to bound 
states. To include the effects of a possible bound state 
we can use the non-perturbative expression (|28|) for the 
scattering amplitudes: given an approximate expression 
for the bound state wavefunction of energy — —Xk^, 
we can calculate the scattering amplitudes and obtain 
corrections from the bound state to the pair and back- 
flow potential in the same way as shown above for the 
scattering states within the Born approximation. In a 
similar way one should proceed to obtain approximations 
for the pair and backflow potentials for systems where 
the interatomic potentials cannot be treated within the 
Born-approximation, for example potentials dominated 
by a hard core. 

Of course, in the case of a single electron in an external 
potential we can solve the Schrodinger equation by other 
means and obtain the "best" pair and backflow poten- 
tials from the exact (numerical) solution. This is done 
below for a perfect crystal using a band structure calcu- 
lation. However, the simple perturbative approach above 
provides an easy way to get some intuition for the pair 
and backflow potential, and is already good enough to 
determine their asymptotic properties. These properties 
are expected to hold in the many-body case: the short- 
range properties are typically determined by two-body 
collisions and the influence of the remaining particles on 
the long-range properties are usually well described by 
an effective single particle potential. Many-body pertur- 
bation theory, which we discuss next, leads to similar 
expressions. 
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B. Many-body perturbation theory 

We now make an expansion of the exact N particle 
wavefunction \(f)) of the interacting system around the 
non-interacting (ground) state |(/)o); the ground state 
without both the electron-electron and electron-proton 
interaction. Let Ofc (a^) be the annihilation (creation) 
operator for an electron of wavevector k. Expanding in 
particle-hole excitations, we have: 



(1 + 



q,ki,k- 



(34) 

The problem is reduced to determining the coefficients 
Q^fci.fc2,i?- Just as in the single particle case, a further ex- 
pansions of aki,k2.,q around fci = or A:2 = together 
with an exponentiation brings the wavefunction into the 
desired functional form, thereby determining the pair and 
backflow potentials. To avoid over counting, we assume 
that the summation in Ea. H34|l goes only over distinct 
states so that it is sufficient to antisymmetrize the wave- 
function at the very end, once we have calculated the 
perturbative corrections. We have limited the expansion 
in Ea. (|34|l to the leading order corrections, particle- hole 
excitations; the generalization to include higher order ex- 
citations is straightforward, but not necessary to calcu- 
late the pair and backflow terms in the wavefunction. 
In order to determine the coefhcients, we write 

I*fci,fc2,g) = a\^^qak^a\^j^^akMo) and multiply these 
states by a constant phase 



afci,fc2., 



kl,k2,q\ 



OC 



ki,k2,q\ 



(35) 



Note that with this phase factor, the right hand side of 
Eq. is given by the expectation value of an oper- 
ator over the true ground state and the coefficients can 
therefore be identified as a TV-particle Green's function^® . 
By considering only particle-hole excitations in Ea. (|34|l . 
the A^-particle Green's function reduces to a connected 
two-particle Green's function and the lowest order mod- 
ifications to the ideal gas ground state of the homoge- 
neous electron gas are therefore related to the two parti- 
cle Green's function G^^^ at equal times^^ or equivalently 



Oikl.k2,q 



-|42^g«fc24i-Kg«fcl 

m) 



(36) 



Summing up particle-hole bubble diagrams (correspond- 
ing to the RPA approximation) results in an effective 
interaction, Uflpyi(p, w). 



ik,Lo) = l-i{p)D{p,Lo) (37) 



where D{p,uj) is the Lindhard function. Perturbation 
theory can now be arranged to be regulaB^S. We note 
that Eq. (|37|l already contains the correct short- and 
long-range limit of the effective interaction. 



function 




r — > oo 




k (X 


V 


e-'/r 








u 




Y 87rnA r 


/ 

V 2nXk'-' 


2Afc2 


y 


yo - y2r2 + 






"k 


V 


2y2 - 


c{rs) 







TABLE I: Asymptotic properties of tlie Jastrow and backflow 
functions for the 3D electron gas. A = /2m, n is tlie electron 
density, y2 ~ COSSr^, and c{rs) ~ 1 + 0.075^^^/(1 + 0.8/?^). 



Neglecting for the moment any contributions from 
plasmon excitations coming from the poles where 
e{kp,LUp{kp)) = 0, we get 

CkiM.q = (1 - "fci-9)"/ci(l - "-fe2+g)"-fe2 (38) 

VRPAjq, £fci - £ki^q) + VRPAjq, £fc2 - gfc2+g) 

X ~ 7 

2(£/ci + £fc2 ~ £fei-g ~ £k2+q) 

where Uk are the occupation numbers of state k in low- 
est order. Expanding around ki = k2 ^ 0, we get the 
Jastrow and the backflow potential. Including the plas- 
mon excitations will give an important long-range con- 
tribution. However, in the simplest approximation, this 
contribution describes only the long wavelength limit cor- 
rectly, and destroys the correct short distance behavior. 
We will circumvent this problem in the next section using 
the formalism of collective coordinates. 

As already shown in the previous section, we expect a 
more general form for the three-body potential. 



) cx det [expikj(ri + Ar^)] e"^^^)" 



-w 



(39) 



with W= 



where 



Wuir) ~ u{r) Wy{r) ~ y{r). 



(40) 



(41) 



For the interactions of the electrons with static pro- 
tons, we can use the static dielectric function e(fc, 0) to 
obtain the effective electron-proton interaction, and use 
directly the results of the single particle perturbation the- 
ory of the previous section with this screened potential. 

The disadvantage of perturbation theory is that one 
gets correct behavior at long and short distances, but it 
does not provide an unique way to interpolate between 
these limits. In Table ^ we summarize the asymptotic 
properties of the pair and backflow potentials for the 3D 
electron gas. 



C. The Bohm-Pines collective coordinate approach 

Instead of replacing the established form for the Jas- 
trow part proposed by GaskellSS by the direct use of Eq. 
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(|38|l . we prefer to improve the RPA form of Gaskell by 
extending it using perturbative formulas. This is most 
easily done within the framework of the collective coor- 
dinate description of Bohm and Pines using additional 
field variables^^. In this approach, the original Hamil- 
tonian of electrons interacting with each other and with 
static protons is extended by an additional boson field 
with generalized momentum variables 11^ coupling to the 
electron and proton density fluctuations 

i k 
k 

where p| {p^) is the Fourier transform of the electron 
(proton) density, — X^i s"**""'' ^.nd and Pk are 
variational parameters. By imposing the extra condi- 
tions Ilfc^' = on the wavefunction, the ground state 
wavefunction of the new extended Hamiltonian will be 
identical to the original one. For a detailed description 
of this approach we refer to the original literatureSi; we 
will only describe the main steps. 

Carrying out the following canonical transformation 

4>old = exp[iS/h](j)new, -5" = ^ {Mkpl + PkPl) Qk, 

k 

(43) 

where Qj. represents the field coordinate conjugate to 11^, 
we obtain an equivalent Hamiltonian 



H = 
where 



-m. 



Hint — 



V ^ 



y2 



m 2m 



-z — MkQkC 



(44) 



(45) 



(46) 



E QlQfc'MfeMfc'k • k'e^C^-'^')"-^- (47) 



(48) 



Now the ground state of the additional field in the 
zero*'' order Hamiltonian, Ea. H44() . is simply given by 
harmonic oscillator ground states of frequencies — 

{nk^M^/my/^, 



Transforming back and applying the subsidiary condi- 
tions replaces the field operator H^ by Mkp% + PkPk and 
the zero*'' order wavefunction is in the Slater- Jastrow 
form 



exp 



1 ^ Mip^kPt + 2MkPkpUpl 



V 



2hVk 



(50) 



up to a constant factor. Instead of using Mk — 
{vk)^^'^0[kc — k) for the long wavelength part up to /cc, 
and optimizing the cut-off fee, as done in the original 



work of Bohm and Pines, we can use 



and 



for the 



electron-electron and electron-proton Jastrow part taken 
in the RPA approximationSi and relate these functions to 
Mfe and Pk- The resulting residual electron-electron and 
electron-proton interaction is screened, since — > Vk 
and MkPk Vk in the long wavelength limit — > 0. 
A second unitary transformation using 



1 



mujp{k){hLjp{0) + ek) 



(51) 



eliminates Hint to first order. Here, ujp{k) is the plasmon 
frequency at wavevector k. Note that this transformation 
brings the wavefunction into the backflow form. Further- 
more, we treat the remaining terms of the Hamiltonian 
perturbatively as shown in the previous subsection. 

The detailed functions we used for the electron gas and 
for metallic hydrogen are given in the appendix and the 
numerical tests are given in Section V. 



IV. COMPARISON WITH THE BAND 
STRUCTURE WAVEFUNCTION 

In this section, we consider another approach of gen- 
erating backflow functions. As in the discussion of the 
single particle perturbation theory in the last section, 
we consider a perfect lattice of protons in which a single 
electron moves. It is straightforward to expand the wave- 
function in plane waves and obtain a precise numerical 
solution of the one electron problem by diagonalization of 
the Hamiltonian matrix. We study to what extent we can 
recast the "band structure" wavefunction into a backflow 
form. The advantage of this approach is that we are eval- 
uating the entire non-linear effect of a lattice of protons 
on the electron wavefunction, or orbital, which for a per- 
fect lattice is a Bloch wave. However, effects of electron 
correlation or screening are absent for this model. 

As was done in Eq. (26), the exact one-electron wave 
function is expanded in plane waves: 



= E 



i(G-l-k)-r 



(52) 



fie^[exp[iki • r^]] exp 



-E 



nl.n. 



V ^ 2hVk 

k 



(49) 



where G is a reciprocal vector of the lattice and k the 
crystal momentum. We then obtain numerical values for 
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Cfc.G by conventional diagonalization of the Hamiltonian 
in this basis. 

First, we study the wavefunction at k = to determine 
the pair part of the wavefunction, U. Neglecting the 
three-body term we have: 

^^.(|r-z,|) = -ln(0o(r)) (53) 

i 

where are the proton positions. Then by fourier trans- 
forming and assuming a Bravais lattice : 

UG = - f dVe-^G"-ln(<^o(r)). (54) 
Jv 

This is shown in Fig. Qand compared to the RPA form 
(solid line) and cumulant form. Note that we only ob- 
tain information about Mk at values of k on the reciprocal 
lattice. It is seen that except for the first few reciprocal 
lattice vectors, the pair wavefunction is determined by 
the cusp behavior. The non-cusp behavior is due to the 
neglect of higher order terms in the cumulant expansion. 
Some effects are picked up by the three-body term of the 
wavefunction. We note that even for the largest lattice 
vector, the values seem to follow a smooth curve, inde- 
pendent of the lattice directions. The k = component, 
though important, will not affect the many-body nodal 
structure or the correlation effects near the fermi surface. 

Now, let us use the same procedure to estimate the 
backflow function. First, we divide out the wavefunction 
at k = to define the backflow functions : 

k- Vn(r) =k-r + iln[</>k/0o]. (55) 

Assuming Yk(r) is the sum of contributions of proton- 
electron terms on a Bravais lattice we get: 

k.Vrfc(r) = l^zG.ky^e'°-. (56) 

G 

Setting these two expressions equal and taking the 
Fourier transform we arrive at: 

Vq^^J ^re^^'"" ln[0k(r)/0o(r)]. (57) 

In general, the function y^' depends on both k and q. 
For small values of k, the ratio approaches a limit, inde- 
pendent of both the magnitude and direction of k. As 
with the pair term, we can only determine ?/k at recip- 
rocal lattice vectors, q and for k in the first Brillouin 
zone. Shown in Fig. (3) is the ratio for several values 
of k evaluated for a simple cubic lattice plotted versus 
q. The dispersion of the values from a smooth curve is a 
test of the extent to which the band structure orbital can 
be cast into the form of a backflow function. Note that 
only for the smallest values of q is the backflow function 
appreciable. At intermediate values of q one does seem 
some effect of " non-backflow" behavior, however it is not 
clear how important these effects are. At large q, we see 
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FIG. 3: The backflow function yq versus the wavevector q in 
atomic units for solid cubic hydrogen lattice at Ts — 1.31. The 
solid symbols are computed using different values of k and 
qin the range of 0.01 to 0.1 using band theory and Eq. (56). 
The solid line is the cumulant approximation: yq — — IGir/g®. 
The dashed line is the backflow function optimized for an 
interacting body hydrogen with a Gaussian form. Dotted 
line is from Eq. (69). 



the behavior yq ~ Sire^ /Xq^ shown as the solid line as 
expected from the results of Section II and III. 

In figure^lis shown the error in the band energies with 
a backflow wavefunction (BF) and the results for hav- 
ing no backflow effects. For the comparison we used a 
BF function yq — yo exp(— 6g) fitted to the low q behav- 
ior. Since yq drops off rapidly w.r.t q, it is primarily 
the effects at small k that are important to describe—. 
By definition the energies are identical at /c = and the 
curvature around k = is exactly put in by the back- 
flow ansatz, at least assuming cubic symmetry. We see 
that the errors in the band energy go as fc'* instead of 
k^ for the non-backflow trial function. However, near 
the band edge there are serious problems because our as- 
sumed form does not have mixing of the bands required 
by lattice periodicity. We expect such an effect to be 
much reduced for a disordered system since such degen- 
eracies will not occur. 

This achieves our goal or showing that the dominant 
band structure effects can be interpreted as backflow cor- 
rections, particularly at small k. This implies that the 
changes in the nodal surfaces due to an external potential 
of protons are well approximated by backflow functions. 
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FIG. 4: The error in the band energy of a single electron in 
a bcc proton lattice for = 1.31 as a function of k (in the 
100 direction) using plane waves (dashed line) and using a BF 
function (solid line). Both approximations are exact at the F 
point since a trial function exact at k was used, but for the 
BF trial function, the error oc fc" while in the PW case (zero 
backflow) the error is oc 



The backflow form is a much more succinct descrip- 
tion of the single body wavefunction than the expansion 
in plane waves. In the introduction, we emphasized that 
this improves performance because we no longer have to 
perform the band structure calculation. However, there 
is also an improvement in speed of calculation of the 
orbitals using backflow. The expansion in plane waves 
can be quite slow, since the accuracy versus number of 
terms decreases quite slowly. In previous work on metal- 
lic hydrogen^, we divided the band structure orbital by 
an electron-proton Jastrow factor as an approximation 
to 4'o{r), and then re-expanded in plane waves. The re- 
sulting expansion is much more quickly convergent in the 
number of plane waves since the cusp at Vep = is in the 
Jastrow factor. It takes the sum of many plane waves to 
recover this non-analytic behavior at = z^. Backflow 
takes this even further by using the fact that near k = 
the wavefunction can be expanded in pair terms with a 
higher-order cusp. These pair terms can be conveniently 
and rapidly computed, since much of the computational 
effort is to map each pair of particles (ee or ep) onto a 
grid value for a table look-up. The distances and grid 
values are then used for all of the pair terms: the po- 
tential, the Jastrow, the backflow and the polarization 



terms. 

The problems concerning degeneracies of the unper- 
turbed plane wave functions near the edge of the Brillouin 
zone are common to all analytical approaches considered 
up to now. Without a separate treatment of (nearly) 
degenerate zeroth order (plane wave) states, neither the 
cumulant method (Section II) nor perturbation theory 
(Section III) are able to produce the resulting energy 
splitting at the band edge. A degenerate case will have 
to be treated by including all of the degenerate states in 
the unperturbed basis. 



V. QUANTUM MONTE CARLO TESTING OF 
TRIAL FUNCTION FORMS 

There are two principal simulation methods used to 
calculate the ground state energies of quantum many- 
body systems: Variational Monte Carlo (VMC) and Dif- 
fusion Monte Carlo (DMC). In VMC, one samples the 
square of the wavefunction, and, in DMC, one uses a 
trial wavefunction and the imaginary-time evolution to 
project onto the ground state. VMC is potentially very 
powerful because one can use any wavefunction, as long 
as one can easily compute its values. One can add cor- 
relation directly to the wavefunction, leading to a very 
compact accurate wavefunction. The resulting integrals 
are similar to that of the classical partition function and 
therefore demand a simulation algorithm to evaluate. 
The disadvantage of the variational approach is that one 
needs to use the right functional space in order to get 
satisfactory properties. Though DMC is much less de- 
pendent on details of the trial wavefunction than VMC, 
however, lacking an exact fermion algorithm, the results 
still depend to some extent on the positions of the node 
(or phase) of the trial wavefunction. 

The most straightforward, and rigorous approach to 
determine the trial function is to propose a definite an- 
alytic form, containing some parameters, a. One then 
uses VMC to evaluate the variational energy Ev{a): an 
upper bound to the exact energy as a function of a. One 
can use various techniques to optimize the parameters 
to obtain the lowest energy, the lowest variance or some 
combination of the two. Variational optimization^^' has 
determined good backflow and three-body trial functions 
for the electron gas in both 2 and 3 D . The disadvantage 
of optimization is that beyond general trends, it is hard 
to extract analytic behavior because of the noisy behav- 
ior of the optimization method and the restriction to a 
limited functional form. 

Here we compare several different trial wavefunctions 
on two systems: the 3D electron gas and metallic hy- 
drogen. We employ three estimators of the quality of 
the wavefunction: the variational energy Ey = (i/jTlip), 
the variational variance = (ipnip) —E^ and the DMC 
(flxed-node) energy. The first two properties are sensitive 
to all aspects of a wavefunction; the variance is particu- 
larly sensitive to short-range structure since the energy 
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wavef unction 




a' 


Ebmc 


1 


SJ 

BF3-0 

BF-A 

BF3-A 


1.0669 (6) 
1.0613 (4) 
1.0611 (2) 
1.0603 (2) 


1.15 (2) 
0.028 (1) 
0.029 (1) 
0.022 (1) 


1.0619 (4) 
1.0601 (2) 
1.0597 (1) 


5 


SJ 

BF3-0 

BF-A 

BF3-A 


-0.15558 (7) 
-0.15735 (5) 
-0.15762 (1) 
-0.15773 (1) 


0.0023(1) 
0.00057 (1) 
0.00061 (1) 
0.00050 (1) 


-0.15734 (3) 
-0.15798 (4) 
-0.15810 (1) 


10 


SJ 

BF3-0 

BF-A 

BF3-A 


-0.10745 (2) 
-0.10835 (2) 
-0.10843 (2) 
-0.10846 (2) 


0.00039 (.5) 
0.00014 (.5) 
0.00017 (1) 
0.00016 (1) 


-0.10849 (2) 
-0.10882(2) 
-0.10888 (1) 


20 


SJ 

BF3-0 

BF-A 

BF3-A 


-0.06333 (1) 
-0.06378 (2) 
-0.06372 (2) 
-0.06358 (1) 


0.000064 (1) 
0.000027 (7) 
0.000045 (2) 
0.000056 (1) 


-0.06388 (1) 
-0.06403 (1) 
-0.06408 (1) 



TABLE II: Energies and variances for the 3D electron gas 
with A'' = 54 unpolarized electrons in Rydbergs/electron. SJ 
means a Slater determinant of plane waves times an opti- 
mized Jastrow factor. BF3-0 are the result of the numerical 
backflow-3body optimization!-. BF-A are the results using 
the RPA Jastrow, Eg. 1591 together with the analytical back- 
flow formula, Eg. 1691 . BF3-A with the additional asymmetric 
3-body wavefunction of Eg. l40l41ll . 



fluctuations are larger. However, the DMC energy is de- 
termined only by the positions of the trial function node, 
not by the "bosonic" part of the trial function. The 
VMC/DMC calculations we performed were standard 
ones^. All calculations are done with Periodic Bound- 
ary Conditions (PBC), equivalent to the F point for a 
band structure calculation in a cubic unit cell. Hence all 
trial functions were real. Though twist-averaged bound- 
ary conditions (TABC)— are useful in reducing size ef- 
fects, tests showed the relative accuracy of various trial 
functions can be determined with PBC using real trial 
functions. 

First, we discuss the results using backflow and three- 
body wavefunctions on the 3D electron gas as shown in 
Table II. The results using analytic trial functions give 
results comparable to the numerically optimized back- 
flow results of Kwon et all™. We find that for < 20 
the analytic wavefunction have a lower VMC energy than 
the numerically optimized wavefunction. This is mainly 
due to the inclusion of the long range part of the back- 
flow potential. For all values of the analytic wavefunc- 
tions have a lower DMC energy, implying a more accurate 
nodal surface than obtained by numerical optimization. 
For Ts = 20 the numerically optimized VMC energy is 
lower than that of the analytic wavefunction, indicating 
that at least the 3-body part of the wavefunction becomes 
inaccurate at strong correlations. 

Now we consider the use of these same trial functions 
for a system composed of electrons and protons. To de- 
termine the properties using the optimization method, 
we used the RPA form for both the electron-electron (ec) 
and electron-proton (ep) u{r). We used optimized Gaus- 



sians for both the backflow and polarization terms: 

77(r) = Aexp(-(r-ro)Vw^)- (58) 

Even though the optimal functions may have a long range 
tail, as shown earlier, the additional energy gained is 
small and we neglect the long-range terms in setting up 
the parameterized trial functions. An additional Gaus- 
sian (with ro — so as not to change the cusp value) was 
added to the pair term. We did not include ee backflow or 
polarization terms in the wavefunction. The resulting 10- 
parameter wavefunction was then optimized to minimize 
a linear combination of its energy and variance. Shown 
in Fig. 121 are the optimized backflow functions compared 
with the cumulant value, with the analytic form and with 
the band structure determination. The magnitude and 
shape are similar, though differences are apparent. 

We compare the results with three other wavefunc- 
tions. The simplest is the SJ-PW functions^, which do 
not contain backflow, three-body terms and the orbitals 
are simple plane waves. We also used optimized Slater- 
Jastrow functions with orbitals from a LDA calculation^. 
Finally shown are various analytic backflow calculations: 
one contains only ep backflow (and 3body), the others 
have alse ee backflow (3body) included. 

Shown in Table III are both VMC and DMC calcu- 
lations of various wavefunctions for metallic bcc hydro- 
gen at Ts = 1.31, a density very close to the molecular- 
metallic transition. While the detailed results depend 
on the number of particles, in general we flnd that 
the SJ-PW function is in error within VMC by about 
15mH/atom while the BF is in error by about 4 mH/atom 
and the LDA trial function by about 2 mH/atom. Within 
DMC the SJ-PW is in error by 6mH/atom and the BF is 
as accurate at the LDA trialfunction within the statistical 
error. This analysis of errors is done with the assumption 
that the LDA-DMC energy is exact. As another indica- 
tion of the quality, the VMC wavefunction variance is 
roughly a factor of 3 smaller with the BF wavefunction 
than with the SJ-PW wavefunction. 

We see that for iV = 16 the DMC backflow results 
are even lower than the LDA function. One reason for 
this could be that = 16 has a degenerate ground state 
for a single Slater determinant; many-body effects break 
the degeneracy. It may be that the current simulations, 
though similar to those of Natoli, broke the degeneracy 
in a more favorable way and thus have a lower energy. 
The A^ = 54 system has a non-degenerate ground state 
at the mean field level, a closed shell, so the results may 
be more typical. Finally degeneracy effects are probably 
less important at A^ = 128 since A^ is larger. 

We have tested the relative importance of including 
ee backfiow in the case of metallic hydrogen. Using ep 
backflow-3body only, the analytical wavefunctions give 
considerably higher energies compared to the numerically 
optimized ones. Including ee backflow in the analytical 
forms, they become comparable. One should note that 
the analytical approaches derive ee and ep backflow at 
the same order of approximation; dropping one of them 
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N 


wavcfunction 


Ev 




Ed mc 


16 


SJ-PW 


-0.4754 (2) 


0.0773(25) 


-0.4857 (1) 




LDA 


-0.4870 (10) 




-0.4890 (5) 




BF3-0 ep 


-0.4857(1) 


0.0317 (5) 


-0.4900 (1) 




BF3-A ep 


-0.4798 (1) 


0.0513 (2) 






BF-A ee+ep 


-0.4850 (1) 


0.0232 (1) 


-0.4905 (1) 




BF3-A ee+ep 


-0.4850 (1) 


0.0227 (1) 






BF-A ee+ep+b 


-0.4878 (1) 


0.0181 (4) 




54 


SJ-PW 


-0.5241(3) 


0.0642 (9) 


-0.5329(1) 




LDA 


-0.5365 (5) 




-0.5390 (5) 




BF3-0 ep 


-0.5331 (6) 


0.033 (1) 


-0.5381 (1) 




BF3-A ep 


-0.5261(1) 


0.0516 (3) 






BF-A ee+ep 


-0.5323 (1) 


0.0222 (2) 


-0.5382(1) 




BF3-A ee+ep 


-0.5325(1) 


0.0214 (1) 






BF-A ee+ep+b 


-0.5353 (2) 


0.0178 (2) 




128 


SJ-PW 


-0.4818 (2) 


0.0656 (23) 


-0.4900 (2) 




LDA 


-0.4962 (2) 




-0.4978 (2) 




BF3-0 ep 


-0.4934 (2) 


0.035 (2) 


-0.4958 (3) 




BF3-A ep 


-0.4846 (3) 


0.059 (1) 






BF-A ee+ep 


-0.4928(2) 


0.030 (1) 


-0.4978 (4) 




BF3-A ee+ep 


-0.4926(2) 


0.029 (1) 






BF-A ee+ep+b 


-0.4947(2) 


0.023(1) 





wavefunction 


Ev 




SJ-PW 
BF3-0-bcc ep 
BF3-0-liq ep 

BF3-0-liq ee+ep 
BF3-A ee+ep 

BF3-A ee+ep+b 


-0.4225(8) 
-0.4418(5) 
-0.4433(8) 
-0.4462(8) 
-0.4430(4) 
-0.4464(6) 


0.0812(4) 
0.0447(7) 
0.0710(10) 
0.0482(8) 
0.0548(2) 
0.052(2) 



TABLE IIL Energies for bcc hydrogen at ra = 1.31. SJ-PW 
means a Slater determinant of plane waves times an optimized 
Jastrow factor. LDA means LDA orbitals times an optimized 
1 body factor and Jastrow factor^, BF3-0 ep means opti- 
mized e-p backflow, e-p polarization and Jastrow. Energies 
are given in hartrees per atom. Periodic boundary conditions 
(F point) and Ewald sums were used, a is the variance per 
electron. BF3-A ep are the analytical wavefunctions using ep 
backfiow-3body only, wheras BF-A ee+ep are results with ee 
and ep backflow; BF3-A ee+ep include also ee and ep 3body 
and backflow, BF-A ee+ep+b uses the same wavefunctions of 
BF-A ee+ep but the electron-proton Jastrow and backflow is 
improved by taking into account the effects of a bound state. 



alone is not justified and might explain the importance of 
including ee and ep backflow in the analytical functions. 
The inclusion of 3-body terms does not noticeably affect 
the energies. This is similar to the results for the electron 
gas at comparable densitiesi^. Since the density is close 
to the transition from metallic to molecular hydrogen 
we tried to improve our wavefunction by considering the 
effects of a simple electron-proton bound state on Jastrow 
and backflow in our analytical formulas (see appendix) 
and found significantly lower energies within VMC. 

Note that when ee backfiow is included, it becomes 
necessary to move all electrons together, and for reason- 
able acceptance ratios one must choose an increasingly 
smaller time step as the system size increases. However, 
the more accurate nodal surface gives both a quantitative 
improvement in properties and a qualitative changes in 
some properties such as Fermi liquid parameters^!. 

We also used the CEIMC^ method to generate a col- 
lection of proton positions appropriate to liquid metal- 
lic hydrogen at 5000K, far above the melting tempera- 
ture of the lattice. Using these configurations we tested 
the accuracy of the same trial functions described above. 
See Table IV. The values marked BF3-0 are obtained 



TABLE IV: Energy and variance of liquid metallic hydrogen 
at rs=1.31, and N = 16. The notation of the trial function is 
described in Table II. The entries marked (bcc) are performed 
with the value of the parameters optimized on the perfect bcc 
lattice. The other entries are optimized over 1000 indepen- 
dent protonic configurations taken at thermal equilibrium at 
5000K. All results are using VMC. 



minimizing local energy and variance for 1000 different 
equilibrium configurations. We compared to the other 
ways of determining the backflow functions, either the 
analytic formulas (see the Appendix) or optimized on 
the lattice. While the optimized BF3 functions have a 
slightly lower energy in some cases, this does not com- 
pensate for the difficultly and reliability of performing 
the optimization. We find that the BF3 wavefunctions 
are about 20mH/atom lower in energy that from the SJ- 
PW at the VMC level, and have a lower variance. This 
comparison shows that disorder weakly affects the deter- 
mined functions at least in this experiment. This sup- 
ports our belief that the BF3 wavefunction is "transfer- 
able" to a variety of protonic configurations. In addition, 
we expect the backflow wavefunction to be more effective 
in the disordered system, since the energy degeneracies 
caused by crystal symmetry of a perfect lattice are not 
present. Comparisons using optimized LDA functions to 
support this hypothesis will be reported in a future pub- 
lication. 



VI. CONCLUSION 

What we have shown in this paper is that ideas from 
perturbation theory can be used to generate an explicit 
trial wavefunction beyond the pair level. This gives us 
both an insight into the form of the many-body wavefunc- 
tion and a more efficient quantum Monte Carlo simula- 
tion for disordered systems. This approach has also given 
intuition on the effect of an external potential on the 
wavefunction, even for a single electron. We have shown 
that one can approximate the band wavefunction (a 3d 
table of numbers for each Bloch wave), with three ID 
functions {u, w, and y) valid for all Bloch waves achiev- 
ing reasonable accuracy. It should be recalled that for 
the electron-proton system, there will be these 3 func- 
tions for the ee interaction and 3 functions for the ep 
interaction. We have found analytical representations of 
these functions accurate throughout most of the phase 
diagram of the electron gas and promising for metallic 
hydrogen. 



13 



An important consideration in Monte Carlo is com- 
putational efRciency. For electron-electron backflow, the 
code runs slower due to having to move all the particles 
together. For electron-proton backflow that is not the 
case. You can still move electrons one at a time since all 
the changes in the Slater matrix are confined to a single 
column; each such matrix value is given by a term of the 
same form as a classical force, allowing it to be quickly 
computed once the ep distance have been computed. Ex- 
pansions of single body orbitals in a plane wave basis can 
be quite time-consuming, especially when pseudopoten- 
tials are not used. 

But the most important advantage of the backflow 
wavefunctions is that the form can be easily extended 
to put in effects of electron-electron correlation on the 
nodes. The outstanding problem in the simulation of 
quantum systems is the "fermion sign problem." If the 
nodal surfaces are accurately approximated, then the 
"fixed-node" method will give accurate results. The 
present work, establishes new analytic properties of the 
backflow functions and thus leads to important progress 
in understanding nodal surfaces. In particular, the ef- 
fect of long-range interactions resulting in perturbations 
to the nodal surfaces is important to establish. Strong 
short-range effects can be captured either by energy min- 
imization or by the nodal release algorithm, which can 
solve for the exact wavefunction for relatively short pro- 
jection times or for small numbers of fermionsi^. Fixing 
the relationship between the long-wavelength collective 
coordinates and the nodal surfaces could be crucial in 
obtaining accurate simulations for fermion systems. 

In the above, we have discussed the use of backflow 
functions for simple metals, using plane waves as the 
reference state. It is straightforward to apply the ap- 
proaches explored here to an insulating state. In that 
case, the reference state will be a determinant of Wannier 
functions, in the simplest case, Gaussians. The backflow 
ideas are applicable for suggesting improvements to the 
resulting Slater- Jastrow function. This will be considered 
in future work. A related problem is how to treat bound 
states in metallic liquid hydrogen in a more accurate way. 

Backflow ideas are also useful at finite temperature. 
In that case we need to know how density matri- 
ces will evolve going from high temperature to low 
temperature^. One knows how to put in backflow at 
high temperature. The challenge is to smoothly interpo- 
late to zero temperature since it is clear that the back- 
flow potential must be a smooth function of temperature. 
In the variational density matrix method^ one uses a 
Hartree-Fock approach with a Gaussian basis to deter- 
mined the evolution of the nodal surface of the many 
body density matrix. The various approaches we have 
described here in particular the Bohm-Pines method, will 
be useful in understanding the temperature dependence. 

Another important problem is to generalize these 
methods to treat electrons with core states. The formal- 
ism should generate good trial functions in the valence 
region and can be used with either all-electron methods. 



or pseudopotentials in that region. We hope that with 
some modification the procedures we have discussed will 
be useful in the core region as well. 
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Appendix: Analytic expressions of the trial 
wavefunction 

In this section we summarize the analytic two-body, 
backflow and polarization functions which describe the 
trial functions. We start from the pair-product (Slater- 
Jastrow) wavefunction based on the RPA approximation, 
using 



and 



1/2 



2nVa 



Eq (1 + 2nVq/eq)^^'^' 



(59) 



(60) 



where £q — Ti^q^ /2m = . Here m is the electron mass 
and n is the electronic density. Using a trial function with 
ee and ep Jastrow factors corresponds to the following 
extended Hamiltonian, Ea. l|44|l . with: 



MqPq 



u'^^u'^q2neq 



(61) 
(62) 



Applying the unitary transformation H51() to the wave- 
function, it generates the backflow potentials. 



2AM2 



ep,int 

yq 



2XMqPq 



t^p('7)('^p(0) 



(63) 
(64) 



where we used ujp{q) = SnXe'^n + 2AkpXsa + e'f, for the 



Fermi vector. 



plasma frequencies and kp is the Fermi vector. The 
screened interaction between electrons, Ea. (|15)l . and be- 
tween electrons and protons, Eq.l^HJ, can be treated by 
perturbation theory. Summing up the particle-hole (bub- 
ble) diagrams, using only the zero*'' order plane waves, 
leads to coefficients ak-i.k2,q for the electron gas, as given 
by Ea. H38() . but with an effective interaction and dielec- 
tric constant: 



iosiq) ^Vq-M, 



e,s{q,uj) = \-vtl,{q)D{q,u) (65) 
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where and D{q,uj) is the real part of the Lindhard func- 
tion. As the pair product form already accounts for plas- 
mons, we do not consider any additional plasmon contri- 
butions to Eq. Expanding Eq. around fci = 0, 
we obtain: 



2g2 eqecff(g,£g) 



(66) 



to obtain this formula we have further approximated 
the sum over occupied (unoccupied) states by [S'((7)]2/4 
where S{q) is the ideal gas structure factor, 



where is the position of the i*'* proton and the sum 
extends over all N protons. As single particle or- 
bital we will take the hydrogen ground state, = 
(7ra^)~^/^ exp(— r/flh), with energy — l/2mal; ^ < 1 
is a normalization taking into account the non-zero over- 
lap between orbitals on different sites. Using Eq. (|28|l we 
obtain for the scattering amplitude 



y/Nna^ p2 (a/ + hTF^ 



(75) 



Siq) = 







3q 




2kF 


\2kp ) 



q < 2kF 
q > 2kF- 



(67) 



The screened electron-proton interaction gives a similar 
term 



ep.sr 



Eqq^ ecff('7, 0) 



(68) 



with the screened electron proton interaction v'^^j{k) = 
Adding these two contributions, the total backflow is: 



Vq 



,,ep _ ,,ep,mt 
Uq Uq 



Vq 



Vq 



(69) 



(70) 



We also performed calculations with an additional ee 
Jastrow function ~ v'^{q) / (eqecs{<l, 0)) but this form did 
not lower the energy. Assuming this form disturbs the 
already correct limiting behavior of the Jastrow part m^'^ 
and it^P for g — > and g — » cxd, we took only the portion 
around the logarithmic singularity at 2kp^ by using the 
following additional Jastrow factor: 



-ee,add 



1 



1 



--ep,add _ 
— 



eoff(g,£g) eoff(0,0) 



1 



ecfr(g,0) eoff(0,0) 



(71) 



(72) 



^ep.add jr^j. ^j-^g total electron- proton Jas- 
trow potential, but only vf^ , since the additional term 



We used u'^ 



;,add 



did not improve the variational energies of the 
electron gas. We used the unsymmetrical form of the po- 
larization with different left and right components given 
by Eq. ^ and Eq. (HIJ: 



w''J'{r) = u''P-'"^\r) w;P{r) = y'^f (r). 



(73) 



Analogous forms were used for the electron-electron part. 

For the case of metallic hydrogen we tried to take 
into account the effects of a possible bound state on the 
electron-proton pair and backflow potential. The sin- 
gle electron wavefunction (j)h considering only one bound 
state can be approximately written by 



(74) 



where we have taken a screened Coulomb interaction 
v{r) = — e^e~'^^^''/r with the Thomas-Fermi wavevector 
k^p — 2fcFe^/7rA, and we have neglected overlap effects 
from different sites. From Eq. H27|) we can finally derive 
the corrections to the pair potential, 



ep,b 



and 



V^(i + (g/2fcTF)2) 

47re^ 

[g2 + {a-^ + kTF)^][eb + ' 



V^(l + (g/2fcTF)2) 
1 

[q^ + (a^^ + kTFf]^[eb + eq 

A 

[q' + {a^' +kTFr][ef, + eq]^ 



(76) 



(77) 



for the backflow potential. We have cut-off the short- 
range part of the corrections by multiplying with [1 -|- 
{q/2kTF)^]~^ in order not to destroy the cusp conditions. 
Since we expect a higher energy for the ground state 
of the screened Coulomb interaction than for the pure 
Coulomb potential, we used ah ~ 2ao and A w 1 in the 
numerical calculations. 

All potentials were split into a short range and long 
range par^^ in such a way as optimize the accuracy for 
a given r-space and k-space cutoff. The short range func- 
tion is evaluated in real space and the long range part 
is then calculated by summing over Fourier components. 
Figure El shows numerical values of rj(r) for the 3D elec- 
tron gas. Comparing with the same figure of Kwon et 
al. where these functions were numerically optimized, 
we see that the short ranged functions are very similar 
for < 10 but different at larger r^. Figure shows 
the three-body contribution of i*'* wavefunction. It is a 
rapidly increasing function of Ts and is somewhat nar- 
rower and more structured than the numerically opti- 
mized form. 



15 



0.4 



0.3 



^ 0.2 

(0 



/ \ 
- I \ 

I \ 

I \ 

-I /'^x \ 



0.1 



\ 

w 



_L 



_L 



1 2 
r/ a 



FIG. 5: The change in the quasiparticle coordinate r'q{r) 
(analytic backflow) caused by an electron a distance r away 
in the 3D electron gas. Graphed is only the short range part 
of rj with A'' = 54. The four figures are for Vs = 1, 5, 10, 20 
from the bottom to the top of the figure. Compare to the 
optimized forms in Fig. (2) in Kwon et ali— 
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FIG. 6: The three-body contribution to the logarithm of the 
wavefunction due to three electrons in the 3D electron gas. 
This is just the short range parts of u'(r) for N = 54. The 
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